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Abstract 

We present a method to evaluate on the lattice the leading isospin breaking effects due to both 
the small mass difference between the up and down quarks and the QED interaction. Our proposal is 
applicable in principle to any QCD+QED gauge invariant hadronic observable which can be computed 
on the lattice. It is based on the expansion of the path-integral in powers of the small parameters (rhd — 
m u )/ Aqcd and a em , where m/ is the renormalized quark mass and a e m the renormalized fine structure 
constant. In this paper we discuss in detail the general strategy of the method and the conventional, 
although arbitrary, separation of QCD from QED isospin breaking corrections. We obtain results for 
the pion mass splitting, - M^ = 1.44(13)(16) x 10 3 MeV 2 , for the Dashen's theorem breaking 

parameter e 1 0.79(18)(18), for the light quark masses, [rh d - m u ](MS,2 GeV) = 2.39(8)(17) MeV, 
[m u /md](MS, 2 GeV) = 0.50(2)(3) and for the flavour symmetry breaking parameters R and Q. We 
also update our previous results for the QCD isospin breaking corrections to the Kit decay rate and 
for the QCD contribution to the neutron-proton mass splitting. 



1 Introduction 



One of the primary goals of lattice QCD is to calculate non-perturbatively hadronic observables at the 
level of accuracy required for phenomenological applications. In the flavour physics sector, for instance, 
the combined efforts of the lattice QCD community resulted in calculations of quantities such as the 
and Kis decay rates with relative overall uncertainties of the order of half a percent (see ref. 1 for a 
recent review). These results have been obtained, in most of the cases, within the isosymmetric theory, 
i.e. by neglecting the difference of the up and down quark masses together with the QED interaction and 
by taking into account the corresponding effects by relying on chiral perturbation theory or on model- 
dependent approximations. At the level of precision presently achieved for some flavour physics observables 
isospin breaking effects cannot be neglected any longer. For example, by neglecting the pion mass difference 
(3%) and the kaon mass difference (1%) a systematic error is unavoidably introduced on the corresponding 
determination of the decay rate or on any dimensional quantity if these masses are used to calibrate 
the lattice. 

In ref. 2 we provided a method to calculate the leading QCD isospin breaking effects, i.e. the ones 
associated with the difference of the up and down quark masses, and we checked the validity of the proposed 
procedure by computing the kaon and nucleon mass difference, the difference of the decay constants ratio 
F K +/F n + with respect to the value of the isosymmetric theory and estimated QCD isospin breaking effects 
on the Kis decay rate. The results just mentioned were obtained by relying on the estimates of QED isospin 
breaking effects, often based on model dependent approximations, provided by other groups. With the 
purpose of removing this approximation we briefly discussed in ref. 2 how order a e rrQ QED corrections 
can be calculated on the lattice. 

In this paper we develop a method to calculate leading isospin breaking effects on the lattice by including 
those associated with QED interactions. These are tiny because very small factors, (rhd — rh u )/ Aqcd 
and a em , multiply sizable matrix elements of hadronic operators. Our approach consists in a combined 
expansion of the lattice path-integral in powers of rhd — rhu and d em . We consider the two expansion 
parameters of the same order of magnitude, (rhd — ^u)/^qcd ~ «em ~ £, and neglect in this work terms 
of 0(e 2 ). In this sense we talk of "leading isospin breaking" (LIB) effects. A great advantage of our 
method with respect to other approaches (see for example refs. [3j|7]) is that, by working at fixed order in 
a perturbative expansion, we are able to factorize the small coefficients and to get relatively large numerical 
signals. For the same reason, we do not need to perform simulations at unphysical values of the electric 
charge, thus avoiding extrapolations of the lattice data with respect to a em . 

The expansion of the lattice path-integral in powers of a em leads to correlators containing the integral 
over the whole space-time lattice volume of two insertions of the quark electromagnetic currents multiplied 
by the lattice photon propagator. These quantities have both infrared and ultraviolet divergences that 
must be removed by providing an infrared safe finite volume definition of the lattice photon propagator 
and by imposing suitable renormalization conditions. In this paper we discuss in detail these issues. 

The main results of the paper are 

M 2 + - M 2 = 1.44(13) (16) x 10 3 MeV 2 , 
[M|+ - M 2 KQ ] QED = 2.26(23)(23) x 10 3 MeV 2 , 
[M 2 K+ - M 2 K o] QCD = -6.16(23)(23) x 10 3 MeV 2 , 

1 Through all the paper we indicate the renormalized couplings with an "hat", for example a e rm to distinguish them from 
the corresponding bare quantities, for example 
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[Ml + -Ml ] 



QED 



[Ml, ~ 



QED 



■Tr-\- 7T U 



= 0.79(18)(18) 



[m d - m u ] (MS, 2 GeV) = 2.39(8)(17) MeV 



rhd 



(MS, 2 GeV) =0.50(2)(3) , 



m u (MS,2 GeV) = 2.40(15)(17) MeV 



m d (MS,2 GeV) = 4.80(15)(17) MeV , 



R(MS, 2 GeV) 
Q(~MS,2 GeV) 



m s - m ud 



m d - m u 



ml — m 2 , 



ml -ml 



(MS ,2 GeV) = 38(2) (3) , 
(MS, 2 GeV) = 23(1)(1) , 



Fr+ I F n + 
. F K /F„ 



1 



QCD 



-0.0040(3) (2) , 



[M n - M P ] QCD = 2.9(6)(2) MeV . 



and have been obtained in the n f = 2 theory. The numbers in the first parentheses correspond to the statis- 
tical errors while those in the second parentheses are the systematic errors, mainly due to chiral, continuum 
and infinite volume extrapolations. The results for the quark masses, F K +/F n + and the neutron-proton 
mass splitting are an update of our previous results obtained for these quantities in ref. 2 . Note that, 
because of the QED interactions, ratios of quark masses of different electric charges are renormalization 
scheme and scale dependent. 

At first order in m d — m u and a em the pion mass difference is neither affected by QCD isospin breaking 
corrections nor by electromagnetic isospin breaking effects coming from the dynamical sea quarks. For 
this reason the result for M 2 + — M 2 is a particularly clean theoretical prediction, though our results were 
obtained by neglecting a quark disconnected contribution to M^o of 0(a em m ud ), see eq. (66). From the 



phenomenological point of view this contribution is expected to be very small, i.e. of the same order of 
magnitude of the other 0(a em [m d — rh u \) contributions neglected in this paper. 

The kaon mass splitting includes both strong and electromagnetic isospin breaking effects. With our 
method these can be conveniently separated by implementing the renormalization prescription discussed 
in detail in section [5] (see also ref. 8 ) and the notation qQ ed ,Q cd a i so means that the corresponding 
numerical result depends upon this prescription. The results for M^+ — M^-o, together with those for 
the Dashen's theorem breaking parameter e 1 and for the light quark masses have been obtained within 
the electro-quenched approximation, i.e. by considering dynamical sea quarks as neutral with respect to 
electromagnetism. 

The paper is organized as follows: in section [2] we discuss the general aspects of our method by assuming 
that the regulated theory retains all the symmetries of the continuum theory. In section [3] we provide an 
infrared safe definition of the lattice photon propagator and discuss a convenient stochastic method to 
calculate electromagnetic corrections to lattice correlators. In section [4] we enter into the details associated 
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with the regularization of the fermion action used in this work. In particular we discuss the issue of the 
determination of the electromagnetic contributions to the critical masses of Wilson quarks. In section [5] 
we discuss all the details needed in order to derive the isospin breaking corrections for a given correlator 
and obtain explicit results for the pion and kaon two-point functions. In section [6] we show our results 
for the pion mass difference. In section [7] we discuss the numerical determination of the electromagnetic 
critical masses of the quarks. In section [8] we discuss the separation of QED from QCD isospin breaking 
corrections to the kaon mass difference and in section [9] we discuss chiral, continuum and infinite volume 



extrapolations of our lattice data. We draw our conclusions in section 10 



2 Electromagnetic corrections to hadronic observables 

In this section we illustrate the general strategy underlying the method that we have devised in order to 
calculate LIB corrections to hadronic observables. The method presented here is a generalization of the 
one presented in ref. 2 and is based on a combined perturbative expansion of the full theorjj^] lattice 
path-integral in the small parameters a em and (rhd — rh u )l ' ^qcd- The two parameters are considered of 
the same order of magnitude, 



m d 



A 



QCD 



0(e), 



(1) 



and contributions of 0(e 2 ) are neglected in the present paper. By using this method it is possible to 
calculate LIB corrections by starting from gauge configurations generated with the isosymmetric QCD 
action. All the details associated with the lattice regularization used in this work will be given in the 
coming sections together with the formulae necessary to compute LIB corrections to specific observables. 

In order to calculate 0(a ern ) corrections to a given physical quantity we have to cope with correlators 
containing two insertions of the electromagnetic current multiplied by the photon propagator and integrated 
over the space-time volume. More precisely, the correction to a given correlator is proportional to 



T{0( Xi )) — ► T / d 4 yd 4 z D^(y - z) (0( Xi )J»(y)r (z)) , (2) 



where T(0(xi)) is the T-product of a certain number of local operators, D /JjU (y — z) is the photon propa- 
gator in a fixed QED gauge and J^(x) is the sum of the electromagnetic currents of all the flavours. There 
are two important issues that have to be addressed in order to give a physical meaning to the previous ex- 
pression. The first, the "infrared" problem, concerns a proper definition of the finite volume lattice photon 
propagator. In section [3] we provide a solution to this problem by discussing in detail how the convolution 
integrals appearing into eq. Q can be calculated numerically. In the remaining part of this section we 
illustrate the "ultraviolet" problem associated with eq. ([2]), i.e. the appearance of divergent contributions 
generated by the contact interactions of the electromagnetic currents. The problem is illustrated here in 
continuum-like notation, i.e. by assuming the existence of a non-perturbative regularization that retains 
all the symmetries of the continuum action (think for example of Overlap lattice Dirac operators), whereas 
we shall enter into the details specific to the lattice regularization used in this paper in the next sections. 
In particular, in section [4] we shall discuss the delicate issue of the cancellation of the linear divergences 
associated with the shift of the quark critical masses induced by electromagnetism. 

2 We call "full" the theory with both QCD and QED interactions switched on and (consequently) with rhd / while we 
call "isosymmetric QCD" or simply "isosymmetric" the theory without electromagnetic interactions and with rhd — triu- 



3 



We are interested in the calculation of the electromagnetic corrections to hadron masses and we do not 
discuss the renormalization of the operators 0(xi). These are needed in order to interpolate the external 



states and, in general, are not QED gauge invariant (see sec. 5.1 for an extended discussion of this point). 
The appearance in eq. Q of ultraviolet divergent contributions associated with the contact interactions 
of the quark electromagnetic currents is understood by considering the short distance expansion of their 
product, which reads 

J^(z)J M (0) ~ ci(x)l + ^c^(x)m / ^ / +c Ps (x)G^G^ + ... . (3) 

/ 

The "counter-term" coefficients ci, and c 9s are divergent quantities that must be fixed by specifying 
appropriate renormalization prescriptions. In particular, the terms proportional to c^, can be reabsorbed 
by a redefinition of each quark mass mj in the full theory with respect to isosymmetric QCD, the term 
proportional to c 9s can be reabsorbed by a redefinition of the strong coupling constant (i.e. of the lattice 
spacing) while the term proportional to c\ corresponds to the vacuum polarization and the associated 
divergence cancels by taking the fully connected part of the right hand side of eq. ([2|. 

In order to take into account the dependence of the parameters of the theory, for example m/, with 
respect to a em and to absorb the divergences originating from electromagnetic interactions, one can include 
in the correlator T(0(xi)) explicit insertions of the corresponding operators, for example ofrffi/jf. To put 
the discussion on a concrete basis, let us consider a generic "physical" observable O in the full theory, 

O(g) = G(e\glm u ,m d ,m s ) = (O) 9 , (4) 
where we have used the following compact vector notation for the bare parameters of the theory 

9 = [e 2 ,g 2 s ,m u ,m d ,m^j (5) 

and where the notation (-) 9 means that the path-integral average is performed in the full theory (see 
section [5J. In the previous expressions we listed the bare mass parameters of the three lightest quarks, 
but the discussion can be easily generalized to include heavier quarks. We have called g s the bare strong 
coupling constant and e the bare electric charge. Note that physical observables are QED and QCD gauge 
invariant and depend on e 2 and g 2 s . Our method consists in expanding any observable 0(g) with respect 
to the isosymmetric QCD result 0(g°) according to 

0(9) = 0(f) + {e 2 ^ + [ 9 l - A + [ m/ _ m 0] JL} (g) 

= {Oy° + AO , (6) 
where 

f=(0,(g o s f,m o ud ,m° ud ,m°). (7) 
The notation (■)^° means that the path-integral average is performed in the isosymmetric theory and the 



9=9 U 



expression of A O in terms of (-) 9 and the appropriate reweighting factor is given in eq. (50). 
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The bare parameters g° of the isosymmetric theory can be fixed independently from the parameters g 
by using an hadronic scheme in order to renormalize isosymmetric QCD, i.e. by performing a "standard" 
QCD simulation, by using a suitable number of hadronic inputs to calibrate the isosymmetric lattice and by 
assuming that isospin breaking effects are negligible. The corrections AO to the physical observables that 
have been used to calibrate the isosymmetric lattice vanish by construction with this prescription while, 
obviously, AO is different from zero for any other predictable quantity. On the other hand, by performing 
simulations of the full theory, the parameters g° can also be fixed by matching the renormalized couplings 
of the two theories at a given scale /i* [8]. More precisely, once the renormalized parameters = 
have been fixed by using an hadronic prescription, the renormalized couplings of the isosymmetric theory 
#i (/x) = Z®({i)g® at the scale /i* are fixed by imposing the following matching conditions 



m°(/i*) = m a (p*) . 



m d {n*) + m u (fj,*) 



(8) 



In this work we rely on this prescription by matching the couplings renormalized in the MS scheme at 
fj* = 2 GeV. 



It is important to realize that a physical observable is a Renormalization Group Invariant (RGI) quan- 



tity, 



(9) 



By using these properties, the perturbative expansion of eq. (|6| can be expressed in terms of the renor- 
malized couplings according to 



om = o a°) + 



de 2 



9s 



d 



m f 



dm f 



O(gi) 



(10) 



From the comparison of the previous equation with eq. Q we find in the differential operator language the 
divergent terms proportional to Z mf /Z^ l and Z g jZ® s that correspond to the short distance expansion 
counter-terms and c 9s respectively. In practice, these counter-terms do appear because the renormal- 
ization constants (the bare parameters) of the full theory are different from the corresponding quantities 
of isosymmetric QCD, the theory in which we perform the numerical simulations. Once the counter-terms 
have been properly tuned, our procedure can be interpreted as the expansion of the full theory in the 
renormalized parameters d em and rhd — rh u . 



3 Non-compact QED on the lattice at 0(a em ) 



In this section we discuss the non-compact formulation of lattice QED, the issues associated with the 
expansion of the quark action with respect to the electric charge and address the "infrared" problem 
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mentioned in section [2j i.e. we provide an infrared safe definition of the finite volume lattice photon 
propagator that can be conveniently used in numerical calculations by working directly in coordinate 
space. 

Non-compact lattice QED has been used also in ref. 13], where the effects of electromagnetism have 
been computed non-perturbatively on the lattice for the first time, and in all the other computations 
subsequently performed (see refs. [4j|7] for recent works on the subject). In practice, the non-compact 
formulation consists in treating the gauge potential A fI (x) in a fixed QED gauge as a dynamical variable. 
The quarks covariant derivatives are then defined by introducing the QED links through exponentiation, 

A^x) — ► E^x) = e~ ieA ^ , (11) 

and by multiplying the QCD links for the appropriate U(l) em factors, 

V+[U,A]i/> f (x) = [E^xW'U^xWfix + ri-Mx). (12) 

In previous expressions ef is the fractional electric charge of the quark of flavour /, i.e. ef is 2/3 for up- 
type quarks and —1/3 for down-type quarks. Given our conventions, exact gauge invariance is obtained 
if the fields are transformed as follows 

t/> f (x) — ► e-^V/Or) , 4>M — ► 4> f {x)e- ie ^ , A^x) — > A^x) + V+A(x) , (13) 

and we define 

V+f(x) = f(x + fi)-f(x), V-f(x)=f(x)-f(x-fi), V M = V ^ V/J • (14) 

We want to treat electromagnetism at fixed order with respect to d em and, to this end, we need to expand 
the quarks action in powers of e. This procedure is performed by starting from the explicit expression of 
the lattice Dirac operator Df[U, A\g\ to be used in numerical simulations and by calculating 

J2Mx){Df[U,A;g\-D f [UA9\}Mx) 

X 

= Y J {{e f e)A,{x)V?{x)+ { ^A,{x)A,{x)T»{x) + ..X , (15) 



where Vj(x) is the conserved vector current corresponding to the quark / while Tj(x) is the "tadpole" 
vertex. Both the conserved vector current and the tadpole vertex depend upon the particular choice made 
for the discretization of the fermion action and we shall provide the explicit expressions for Vf and 
corresponding to the regularization used in this paper in the following sections, see eqs. (37). Note that 
tadpole insertions, a feature of lattice discretization, cannot be neglected because these play a crucial role 
in order to preserve gauge invariance at order e 2 . The electromagnetic current and the tadpole vertex to 
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be inserted in correlators are the sums over all the quarks of and with the corresponding charge 
factors, 

f f 

t»\x) = j2M 2jl f( x ) = E( e / e ) 2 ^ r T^]^w- ( 16 ) 



From the validity on the lattice of exact gauge Ward-Takahashi identities (WTI) and from the fact 
that the fermion action is by construction renormalization group invariant, it follows that the expansion 



of eq. ( 15 ) is perfectly well defined and that it can be re-expressed in terms of the renormalized electric 
charge and gauge potential fields by replacing e —> e/Z e and — > A^Z e . Furthermore, at the 0(a em ) at 
which we are working, there is no need to renormalize the electric charge, a problem that has to be faced 
instead at higher orders. 

Once the fermion action has been expanded, the leading QED corrections to a given lattice correlator 
are obtained by considering the time product of the original operators with two integrated insertions of the 
combination A fI (x)J^(x) or with a single integrated insertion of ^2 A fI (x)A fl (x)T^(x). As anticipated in 
section [2] the corrected correlator is expressed in terms of the photon propagator. To give an example, let 
us consider the electromagnetic corrections to the kaon two-point correlator. Among other contributions 
discussed in detail in the following sections, in this case one has to calculate 



= es6 " e2 ( T,M*)Mv) T<0| [«76*](t) V?{x) V:(y) [s 75 u}(0) \0)J 

= e s e u e 2 Y, D ^ x ~ v) T <°l l n ^s](t) V»{x) V^y) [s l5 u](0) |0) , (17) 



where the notation (-) A represents the path integral average over the gauge potential A^ (see eq. (48)), 
D^ u (x — y) is the lattice photon propagator and we have ignored the quark disconnected contributions 
coming from the contractions of the vector currents Vf(x) and V^(y) among themselves. 

In order to define the lattice photon propagator we start by considering the lattice action of the QED 
gauge field in Feynman gauge 

SgaugelA] = \ £ A^x) [-VJV+] A^x) = \ A X k ) [^(K/2)] 2 A^k) , (18) 

x,n,v k,/j,,v 

where A^x) is a real field while A^ (k) denotes its Fourier transform that is a complex field satisfying 
the condition A*(fc) = A^—k). In the previous expression we have explicitly shown the QED action in 
momentum space to highlight a well known problem with the definition of the lattice photon propagator, 
i.e. the infrared divergence associated with the zero momentum mode. The A^ propagator is defined as 
the inverse of the kinetic term and, in order to define the inverse of the lattice Laplace operator — V~V+, 
one has to provide a prescription to cope with its kernel. Any "derivative" gauge fixing condition does not 
constrain the zero momentum mode of the electromagnetic gauge potential, 

V-[A^x) + c] = V-A tl {x) , (19) 
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and, as a consequence, the gauge fixing has to be "completed" by giving a prescription to regularize the 
associated infrared divergence. 

One possibility, widely used in the literature after the original proposal made in ref. 3 , is to make the 
zero momentum mode to vanish identically by sampling the gauge potential in momentum space. It can be 
shown that this prescription results into finite volume effect on physical observables, see section [9j We also 
follow this strategy and set A(k = 0) = with the difference that we work directly in coordinate space, 
thus avoiding Fourier transforms, and calculate the infrared regularized photon propagator stochastically. 
More precisely, by introducing the operator P x projecting a given field on the subspace orthogonal to the 
zero momentum mode, 



(20) 



we calculate the regularized photon propagator 



-VpV+ 



(x - y) , 



(21) 



by following the procedure outlined here below: 



we extract four independent real fields B fI (x) distributed according to a real Z2 noise, 



(B^(x)B u (y)) = 5^ S(x - y) ; 



(22) 



for each field B^x) we solve numerically the equation of motion in Feynman gauge, 



[-VrV+]C M [ J B;x]=P ± B^x); 



(23) 



the solution is 



C^B-x] = 



-VpV+ 



B v {x) = Y,D^{x - z)B v {z) 



(24) 



and the field C^[B] x] is a functional of B^ 
• the photon propagator is thus obtained by using the properties of the Z2 noise according to 



{B„(y)C„[B;x}) B =J2 D » P ( X ~ z ) {B»{y)B p {z)) B = D^(x - y) 



(25) 
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This procedure relies on the actual possibility of obtaining a numerical solution of eq. (23). By working 
in double precision we have been able to obtain a stable and efficient numerical bi-conjugate gradient 
stabilized (bicgstab) inverter. The solution is obtained with about hundred iterations on lattice volumes 
as large as V = 96 x 48 3 . Coming back to the example discussed previously, we have that the infrared 



regularized version of eq. (17) can be calculated on the lattice according to 
- <jt> = E D U X - V) T (°\ [«76*](*) V?{x) V:(y) [s l5 u}(0) |0> 

x,y 

= (Y,Mx)CAB;y] T(0\ [u l5 s}(t)V^(x)V:(y) [s l5 u}(0) |0> 



= -(j^B^C^y] Tr{~f 5 S s [U;t - x\T^S s [U;x\ lh S ud [U;-y\T v v S ud [U;y -t]}\ , 

\ x,y I 

(26) 

where we used the compact notation Sf[U] to indicate the isosymmetric QCD lattice quark propagator 
Sf[U;g°] obtained by inverting the Dirac operator Df[U] = Df[U,g°] (see eq. (31) below). The problem 
of the numerical calculation of the diagram appearing on the left hand side of the previous expression is 
thus reduced to the calculation of two sequential propagators. More precisely, one can solve the following 
two systems 

{D f [U]¥ B }(x) = Y, B ^ x )KSfP;x}, 

{D f [U]¥ c }(x) = J2 C ^ x ^vSf[U;x], (27) 



for different values of the B^(x) and C^[B\ x) fields (we have used 3 electromagnetic stochastic sources per 
QCD gauge configuration) and then calculate the corrected correlator according to 



606716 



<^> = -e s e u e 2 < Tr { [^(t) ** B (t) } f . (28) 



In the previous expressions we have been assuming that the lattice Dirac operator, and consequently the 
conserved vector current, satisfies the property Sf[U ; x]^ = j^Sf[U ; —^75. Generalizations of the previous 
procedure to calculate all the other correlators (fermionic Wick contractions) appearing in this paper can 
be readily obtained. 



4 Fermionic lattice action 

In this section we enter into the details of the fermionic lattice action used in this work, namely the 
maximally twisted Wilson action 9, 10 . In order to minimize cutoff effects and statistical errors we have 
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been working within a mixed-action approach |11||12| . In particular, the results described in the following 
sections have been obtained with the action S = S sea + S va i. The sea quark action is given by 

Ssea = ^2{q u D+[U,A}q u + q d D-[U,A]q d } , (29) 

X 

where qf are fermionic variables and the [U, A] lattice Dirac operators are 

Df[U,A] t/>(x) =m f iP(x) ± i l5 (mf + 4)t/>(x) - ^ ±% ^~ 1r Ufl (x)[E^x)] e fxP(x + M ) 

- E ±f V 1 " u ^ x - Hi* - - • 



2 



(30) 



Note that the operators [?7, A] depend upon the bare parameters of the full theory and we used the 
- - ^ • ^ ^/[^j ^] = ^/[^j ^5 The corresponding operators [[/] = D^[U ; #°] of the isosym- 



metric theory are 



Df[U}^(x)=m° f ^(x) ± i 75 K r + 4)^) - E ^ 75 2 — + 



(31) 
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Concerning the content of the valence sector, we have considered a doublet of fermionic fields for each 
flavour, ipj = (il)~f ,i/jJ), and a corresponding doublet of bosonic fields (pseudo-quarks), 0j = (0jr,0~). 
The fields within the same doublet have the same mass ra/, the same electric charge ef but opposite 
chirally rotated Wilson terms. Calling r % the Pauli matrices acting on a flavour doublet and defining 

D f [U,A}= 1 -^D+[U,A]+ 1 -^Dj[U,A], (32) 

we can write the action for the matter fields in the valence sector in the following compact notation 

S va i = ^{^fDfp.A^f + ^fDfp.A^f} . (33) 
i> 

As far as the mass splitting of the pions (or of the nucleons) is concerned, the mixed action setup used 
in this paper allows to compute observables with 0(a 2 ) cutoff effects at the price of introducing unitarity 
violations that disappear when the continuum limit is performed (at matched sea and valence renormalized 
quark masses the resulting continuum theory is unitary). For each correlator, by possibly replicating some 
of the valence matter fields, the choice made for the action allows to consider only the fermionic Wick 
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contractions that would arise in the continuum theory, thus avoiding the introduction of (finite) isospin 
breaking lattice artifacts. The resulting diagrams are then discretized by using for each quark propagator a 
convenient choice of the sign of the twisted Wilson term. In practice we consider for any meson interpolating 
operators of the form 

H =iP+r<Pj 2 . (34) 

The resulting correlators have reduced cutoff effects and smaller statistical errors with respect to the other 
possible choices of 0#, as for example t^T'*/^, see refs. |IlJ[l2j. In the case of the connected fermionic 
Wick contraction arising in the neutral pion two-point functions we use O^° nn = (u+j 5 u~ — d + ^ b d~) / \/2. 

In the following we have also computed isospin breaking corrections to the kaon masses. In this case the 
strange quark is quenched and the theory violates unitarity also in the continuum limit. In the calculation 
of the kaon mass splitting additional violations of unitarity will be introduced when we shall consider what 
we can call the "electro-quenched" approximation. This approximation consists in forcing the sea quarks 
to be neutral with respect to electromagnetic interactions and is implemented by replacing S sea with 

StTa = Y,{^ D ul U ]Qu + q d Da[U]q d } . (35) 



4.1 Critical mass counter— terms 

We now discuss the problem of the (re) tuning of the critical masses necessary in the presence of electromag- 
netic interactions when the lattice fermionic action includes a Wilson term. In this case eq. ([3| is modified 
both on the left-hand side, to take into account the presence of the tadpole vertices of the different quarks, 
and on the right-hand side, because of the appearance of additional divergent contributions that have to 
be reabsorbed by a redefinition of the critical masses. The lattice version of eq. Q corresponding to the 
regularization used in this paper is thus given 

J"(a;)J M (0) + 

~ ci(ar)l + J2 c k( x ^f i ^r 3 ^f + ^ c^m/V^/ + c gs (x)G^G^ + 
/ / 

(36) 

where (see eqs. ([l6|) above) the exactly conserved vector current and the tadpole vertex corresponding to 



the Dirac operators of eqs. (30) and (32) can be obtained by expanding the lattice quark action according 



to eq. ( 15 ) and are given by 



3 The twisted lattice action, even in the presence of electromagnetic interactions, enjoys the symmetry P x T>d x (mj i— Y 
-TUf) where P is the ordinary lattice parity and T)^ is the lattice version of the transformation that replaces a generic operator 



O(x) of dimension d with e ldlT 0(— x). It follows (see ref. |10||13] ) t hat parity-even physical observables are automatically 
O(a) improved and that the operator G^ V G^ V does not appear in eq. ( |36| >, though cutoff effects proportional to the insertions 
of (certain combinations of) the parity-odd operators rhfG^ v G^ v are not forbidden by the symmetries of the theory. 
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T?(x) = $ f (x) lT \ 7jJ U^x^fix + $ } {x + fi) lT \ + ^ Ul(x)^ f (x) . (37) 



Note the presence in eq. (36) of the critical mass counter-term coefficients c£(x). In case of standard 
(untwisted) Wilson fermions one would get a similar expression with the scalar operators ipfi/jf instead of 
the pseudoscalar parity-odd operators ipfij^T 3 ipf. 

We have considered two different strategies to determine the counter-terms associated with the elec- 
tromagnetic shift of the critical masses, both based on the use of the WTI of the continuum theory. The 
first strategy can be used with both standard and twisted Wilson fermions and is based on the Dashen's 
theorem, a consequence of chiral WTI valid in the massless theory, i.e. the theory with 

rhf {m n , m d , m s } . (38) 

According to the theorem, even in the presence of electromagnetic interactions, the neutral pion and the 
neutral kaons are non-singlet Goldstone's bosons, so that 

m/ = M n o=M K o=0. (39) 

Furthermore, from the vector flavor symmetries which remain valid in the massless theory even in the 
presence of electromagnetic interactions it follows that 

m f = M 7r+ =M K+ (40) 

and that the critical mass parameters of the down and of the strange quarks are equal. By using this 
observation, after a detailed discussion of the corrections to kaon and pion masses, we shall give explicit 
formulae to determine the critical mass counter-terms by imposing the validity of eqs. (39), see section [7] 



The second approach to tune the critical mass parameters is peculiar to chirally twisted Wilson fermions 
and is commonly used to implement the maximal twist condition in simulations of isosymmetric QCD, see 



ref. 14 . By starting from the explicit expression of the lattice Dirac operator for a given valence flavour 



doublet, eqs. (30) and (32), one realizes that the critical mass of each valence quark can be separately 
tuned by working in the theory with electromagnetic interactions and with massive non-degenerate quarks 
(m s > rhd > m u > 0) and by imposing the validity of the following vector WTQ 

W f (g) = V M < [trfr^f] Or) [VvtV^] (0) >* = , / = {u, d, s} . (41) 

Also in the case of this WTI we shall derive explicit formulae corresponding to its expansion in powers of 
e in section [3 



4 A similar procedure to tune the critical masses can be used with standard Wilson fermions by starting from the PCAC 
relation and by using the knowledge of the PCAC quark mass in the isosymmetric theory. 
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Before closing this section we discuss the generalization of eq. (|6| required to take into account the 
dependence of a generic lattice observable on the critical masses m^ r . More precisely, by enlarging the 
parameter space of the full theory, 

0(g) = 0(e 2 ,g^m u ,m d ,m s ,m c u r ,m c d r ,m c s r ), g = (e 2 , g 2 s ,m u ,m d ,m s ,m c u r \mj ,raf) , (42) 

and by calling mg r the single critical mass parameter of the symmetric theory, we see that isosymmetric 
QCD simulations correspond to 



9° 



(0, (g° s ) 2 , m° ud , m° ud , ml m% , m c J , m c j) . (43) 



The value of mg r has been precisely determined in re f. |14| in the isosymmetric theory by requiring the 
validity of the vector Ward-Takahashi identity of eq. Ell) with m® = 

W ud (f)=0 — ► mZ . (44) 

Our gauge ensembles have been generated at this well defined value of critical mass for each f3° = 6/(g!?) 2 
(see Appendix |A| and, for this reason, we cannot tune the value of the different mj by looking at the 
dependence of TV/, or of M n o and M^o, on the quark critical masses used in the simulations. On the 
other hand, the LIB corrections to any observable can be obtained by making an expansion, at fixed 
lattice spacing, with respect to the differences mj' — mg 7 " which represents a regularization specific isospin 
breaking effect induced by the electromagnetic interactions. The generalization of eq. (|6| to be used on 
the lattice with Wilson fermions is 



a ° - { e2 ^ + w-«» 2 ]4 +K - m?i ^ +im ' r - msri ^f} o<a 



(45) 



In the next section we discuss in detail how eq. (45) can be used to expand the lattice path-integral and 



to derive explicit formulae for the calculation of the LIB corrections. 



5 Expansion of the lattice path-integral at 0(a em ) 

In this section, by following the strategy outlined in the previous sections, we discuss the details concerning 
the derivation of the formulae necessary to calculate the LIB corrections to specific observables. The 
starting point is the path-integral representation of the observable in the full theory, 

fdAe- s ^W dUe-P s '™'*M Y\ n / det (Df[U, A; g\) 0[U,A;g] 

O(g) = (Oy = , , (46) 

JdAe-SwM dU e-P s 9^ 9e [U] det (Df[U,A;g\) 



where S gauge [A] has been given in eq. (18) and is a functional of the gauge potential A^ S gauge [U] is the 



QCD gauge action (f3 = 6/g 2 ) and is a functional of the link variables U^x), D f [U,A;g\ are the Dirac 
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operators defined in eq. ( 30 ) . Note that when the masses and the charges of the light quarks are different 
the product of the determinants of the up and of the down is not positive definite unless one resorts to 
lattice regularizations in which the determinant of each quark is separately positive. We want to express 
the observable O(g) in terms of the path-integral average in the isosymmetric theory, i.e. 



o{f) = {oy = 



JdUe-P° s ^M l\y =1 det (pf[U;?]) 0[U] 
JdU e-0 os »<™««M ri/iidet (Df[U;g°]) 



(47) 



where the Dirac operators D^[U;g°] have been defined in eq. (31). This can be done by introducing the 



appropriate reweighting factor and the functional average (•) with respect to the free photon field, 



R[U,A;g\ = e 



r[U,A;g\ 



n f n f det (Df[U,A;g\) 

iU,A;g\ = l[r f [U,A;g\ = T] V^T 

/=i 7=i det (pf[U;g°]) 



( )A ^ JdAe s — 



[A] 



0[A] 



JdAe- 



•-[A] 



(48) 



Eq. ( 46 ) can be conveniently rewritten as follows 



(RO) 



A, if 



(R) 



A,g° 



R[U,A;g\ 0[U,A;g\ ) J 



(R[U,A;g\ ) 



A \9- 



(49) 



and leading order isospin breaking corrections can now be obtained by applying the differential operator 



A defined in eq. ( 45 ) to the observable O defined in eq. ( 49 ) . More precisely 



AO = {A{RO)) A ' r - (AR) A ' f (of = 
(AO[U, A;g\\ s=go ) A ' f + {(A(RO - O) [U,A;g\\ s=g0 ) A '' g 



{AR[U,A;g\\ g=s0 ) A ^(O[U;g°}f}. 



(50) 



In the previous expression we have separated the term (AO) A ^° , representing the correction to the given 
observable, from the contributions in curly brackets coming from the corrections to the reweighting factor 
and, consequently, to the sea quark determinants. In the following we shall call these contributions "vacuum 
polarization terms" or "disconnected terms" . 
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Once the quark fields have been integrated out from the path-integral, as implicitly done in the expres- 
sions above, in order to calculate the LIB corrections we must be able to apply the differential operator A 
to the Dirac operator and to the quark propagator. To this end, it is useful to observe that 



d{0) (e 2 ) 



d{e 2 ) 



; 2 =0 



ld 2 0[A;, 



de 2 



e=0 



(51) 



and to consider the following expressions and related graphical representations 



l d 2 S f 
2 de 2 



dD f dD f 
Sf^Sf—Sf 



de 



1 tPDj_ 

2 f de 2 



Sf 



9S f 
drrif 



dm c . 



= -S 



^ dm f ^ 

± dD t s± 

f dm c r r f 



(52) 



The graphical representation given in the last of the previous formulae, corresponding to the derivative of 
the quark propagator with respect to the critical mass, is specific to the lattice Dirac operators used in 
this work and the =p signs correspond respectively to defined into eq. (30). In the case of standard 
Wilson fermions red and grey "blobs" would coincide. All the disconnected contributions coming from the 



reweighting factor can be readily obtained by using eqs. (52). For example, 



OR 



Q 



[U} = 



2 de 2 



2 V de 



TV S 



dD t 



V 



de 



(53) 



In writing eqs. (52) and (53) we assumed that the derivatives have been evaluated at g = g° and that the 
functional integral (-} A with respect to the photon field has already been performed. Note however that, 



in order to apply the operator A to the product (R[U,A;g\ 0[U,A;g\) (see eqs. (50) and (51) above), at 
fixed QED gauge background one also needs the following expressions for the first order derivatives of the 
quark propagators and of the quark determinants with respect to e 



de ~ bf de bf ~ €f 



drj_ 
de 



= Tr 



(54) 
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A concrete example of application of the formulae given in eqs. (52) and (53) is represented by the 



correction to the Sy quark propagators worked out below 



(e/e) 2 



(e/e) 2 



[rrif — m°f] (x) =p [m 



2^ e h — * — + e Z^ e h e f2 — * — 



£ ±[m£ - mgn > + J2 \m h ~ m° fl ] > + [ 9 l ~ (rf) 2 ] 



(55) 



Here quarks propagators of different flavours have been drawn with different colors and different lines. 

The formulae above have been explicitly displayed not only because they represent the building blocks of 
the derivation of the LIB corrections to the hadron masses discussed in the following, but also for illustrating 
the implications of the electro-quenched approximation (see eq. (35) above). This approximation is not 
required in the calculation of the pion mass splitting because the quark disconnected diagrams containing 
sea quark loops are exactly canceled in the difference of AM 7r + and AM^o (see eq. ( [66] ) below). This does 
not happen in the case of the kaon mass difference, see eq. (69). Quark disconnected diagrams are noisy 
and difficult to calculate and, for this reason, we have derived the numerical results for M K + —M K o within 
the electro-quenched approximation. The perturbative expansion of the electro-quenched theory, i.e. the 
theory corresponding to the action Sf^ for the sea quarks, is obtained in practice by setting g s = g° s and 



r f [U,AJ ] = l. 



(56) 



In the electro-quenched approximation all quark disconnected contributions are absent. It follows that in 



this theory eq. ( 55 ) simply becomes 



(e/e) 2 



[m f 



T[mf 



(57) 



5.1 LIB corrections to hadron correlators 



In order to extract the mass of a given hadron H, by including electromagnetic interactions and QCD 
isospin breaking corrections, we start by considering in the full theory the two-point correlator of an 
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interpolating operator On(t^p = 0) having the appropriate quantum numbers, 



C HH (t;g) = < H (t) OUO) ) S = Z H e~ tM t 

e M H _ ChhH - l;g) _ 
CHH{t;g) 



(58) 



where the dots represent non leading exponential contributions to the correlator. It is important to stress 
that, if H is an electrically charged particle, the correlator Chh{^ g) is not invariant under U(l) em gauge 
transformations. For this reason it is not possible, in general, to extract physical informations directly from 
Zh, the residue of the pole corresponding to the hadron H (see also ref. [l6] concerning this point). On the 
other hand, the mass of the hadron Mh is gauge invariant and finite in the continuum limit, provided the 
parameters g of the action have been properly tuned. It follows that, at large times and at any given order 
in a perturbative expansion in any of the parameters of the action, the ratio Chh{P — 1] g) / C h h (t m , g) is 
both gauge and renormalization group invariant (up to discretization effects and exponentially suppressed 
contributions). From eqs. (58) it follows 



C H H(t;g) = C H H(t;f) 



AC gg (t) 
C HH (t;g°) 



AM H = M, 



H 



= -d t 



AC gg (t) 
C H H(t;g°) 



(59) 



where we have defined 



d t f(t) = f(t)-f(t-i) 



(60) 



and A is defined in eq. (45) 



In our lattice simulations we have enforced periodic (anti-periodic) boundary conditions for the gauge 
(matter) fields along the time direction. For this reason, we have extracted the correction to pseu- 
doscalar meson masses by fitting the ratio ACpp(t)/Cpp(t; g°) of corrected over uncorrected pseudoscalar- 
pseudoscalar two-point functions according to the following functional form 



ACppjt) 



const. + AM P (T/2 - t) tanh [M P (T/2 - t) 



(61) 



where the constant term contains the correction to the residue of the pole corresponding to the lightest 
state of mass Mp and T is the extension of the time direction of the lattice. The formula above is obtained 
by noting that a pseudoscalar-pseudoscalar correlator is even under the symmetry t \-> T — t and by using 
ordinary perturbation theory in order to predict the time dependence of corrected correlators, see ref. [2] 
for further details concerning this point. In the following we continue to use the symbol dt but, when 
referred to lattice correlators, we actually mean the operation that allows extracting the coefficient AMp 



by a fit of the numerical correlators according to eq. (61). 
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As explained in section [4j in order to minimize cutoff effects and optimize the numerical signal, we work 
in a mixed action setup and extract both charged and neutral meson masses from two-point correlators of 
twisted Wilson quarks having opposite chirally rotated Wilson terms. In practice, the use of interpolating 
pseudoscalar operators of the form ^~f 1 7 6fl Pj 2 ( see ec l- (^4) above) corresponds to 



+ 



(62) 



The same assignment of Wilson parameter signs has been used also when the two quark propagators 
correspond to the same physical flavour (for example = f 2 = u), see appendix [A] for further details. 

5.2 Pion two— point functions 

Given the observations made in the previous subsection, we now derive the leading isospin breaking cor- 
rections to pion masses by using the same technique employed to obtain the corrections to the quark 
propagator. In the case of the charged pions we can start from the full theory correlator 



C w+w -(t;g) = ( [u+ 1 b d-]{t,p = 0) [d-j 5 u + }(0) > 



(63) 



and apply the differential operator A defined in eqs. (45). We get 



AM. 



e u e d e 2 d t — (e 2 u + e 2 d )e 2 d t 





2[m ud - m° ud ]d t 



(e u + e d )e 2 ^ e f d t 




f=sea 



(m^ r + m c J — 2rriQ r )dt h [isosym. vac. pol.] 



(64) 



where m u d = (rrid + m u )/2 is the bare isosymmetric light quark mass. In the case of the neutral pion we 
obtain 



AM T o = - 



el + e 2 



-e 2 d t 



<s> 



(el + e\)e^ 





2[m ud - m° ud )d t 



18 



(e u + e d )e 2 ^ e f d t 




f=sea 



(mZ + mJ-2m^)d t 



(e M - e d f 



e 2 d t 



66 



[isosym. vac. pol.] 



(65) 



The sea quark propagators have been drawn in blue (and with a different line) and the isosymmetric 
vacuum polarization diagrams have not been displayed explicitly. By combining the previous expressions 
we find the elegant formula 



(e n - e d ) 2 



<S>-C 6 



(66) 



All the isosymmetric vacuum polarization diagrams cancel by taking the difference of AM 7r + and AM^ 



together with the disconnected sea quark loop contributions explicitly shown in eqs. (64) and (65). Note 



in particular, the cancellation of the corrections/counter-terms corresponding to the variation of the sym- 
metric up-down quark mass m ud — m^ d and to the variation of the strong coupling constant g 2 — (g^) 2 . 
This is a general feature: at first order of the perturbative expansion in a em and rhd — rh Ul the isosym- 
metric corrections coming from the variation of the stong gauge coupling (the lattice spacing) , of m u d and 
of the heavier quark masses do not contribute to observables that vanish in the isosymmetric theory, like 
the mass splitting M 7r + — M^o. Furthermore, as already stressed, the electric charge does not need to 
be renormalized at this order and, for all these reasons, the expression for the pion mass splitting can be 
considered a "clean" theoretical prediction. 



On the other hand, the lattice calculation of the disconnected diagram present in eq. (66) is a highly 



non trivial numerical problem and we shall neglect this contribution in this paper. Relying on the same 



arguments that lead to the derivation of the flavor SU(3) version of the Dashen's theorem, see eq. (39), it 
can be shown that the neutral pion mass has to vanish in the limit rh u = rhd — f° r arbitrary values of 
e n , ed as well as the masses rhf and the electric charges ef of the heavier quarks. This happens because 
the electric charge operator is diagonal in the up-down space and commutes with the isospin generator r 3 . 
Once the critical mass counter-terms m c ^ d — mg r have been properly tuned, the contributions to eq. (65) 



can be separated by the dependence with respect to e n , and e/ of the different coefficients. It follows 
that the disconnected diagram of eq. (66) vanishes in the SU(2) chiral limit and, consequently, it is of 
0(a em m u d) . Neglecting this 0(a em m u d) diagram we are thus introducing a small systematic error that, 
from the phenomenological point of view, can be considered of the same order of magnitude of the other 
0(a em [rhd — rh u \) contributions neglected in this paper. 
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5.3 Kaon two— point functions 



By repeating the analysis performed for the pions in the case of kaon two-point functions, we obtain the 
following result for the LIB corrections to M K + 



AM K + = 




+ [isosymmetric vac. pol.] , (67) 

where the strange quark propagator is drawn in red and the sea quark propagators in blue (with different 
lines). The leading corrections to M K o are given by 



AM K o = 




+ [isosymmetric vac. pol.] . (68) 
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By taking the difference of the last two expressions we get 



M K+ -M K o = (el-e 2 d )e 2 d t -^—-(e 2 u -e 2 d )e 2 dt 




o 

O O Gb 

O O ' o 

(69) 



where we defined 



Am ud = md 2 W " , Amf = mf - m% , (70) 

and used the relation e s = — (e u + e^). Also in the kaon sector, by taking the difference (AM K + — 
AM^o), all the isosymmetric vacuum polarization diagrams cancel as well as the corrections/counter- 
terms corresponding to the variation of the symmetric up-down quark mass m u d — m^ d and of the strange 
quark mass m s — m®. The "sea-tadpole" diagrams don't cancel and contribute to M K + — M K o. These 
terms, absent in the case of the pion mass difference, vanish however in the SU (3) chiral limit and/or within 
the electro-quenched approximation that we shall employ in sections [7] and [8] to obtain our numerical results 
for the kaon mass difference. 



6 Pion masses 



In this section we discuss our results for the physical pion mass splitting. The starting point is eq. (66) 
that, by neglecting the disconnected diagram coming from the neutral pion, can be conveniently rewritten 
as 

K xch (t)= — > M 2 + -M 2 = (e u -e d ) 2 e 2 M 7T d t Rl xch (t) . (71) 



In the left panel of Figure [T] we show the mass of the pions in the isosymmetric theory, M^, as extracted 
from the unperturbed correlators C^^^g ). The data are shown for different values of the lattice spacing 
and for different values of the symmetric light quark mass m^ . In the right panel of Figure [l] we show 
the fits of the ratio of correlators R^ xch (t) according to eq. (61). As can be seen, we are able to obtain a 
good numerical signal by using three photon stochastic sources per QCD gauge configuration. The data 
follow the expected behavior as a function of the time variable. 
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t/a t/a 



Figure 1: Left panel: extraction of the pion mass in the isosymmetric theory from Ctttt^; g°). Right panel: fits of R^ xch (t) 
according to eq. ( |6l| . The dark magenta points correspond to (3 = 3.90 and (amy) = 0.0030, the green points to (3 = 4.05 
and (am U( i) = 0.0060 while the blue points correspond to (3 = 4.20 and (am U( i) = 0.0065 (see Appendix |a|. Data are in 
lattice units. 
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Figure 2: — M^ as a function of rh u d as extracted from lattice correlators and converted in physical units by using 

a . Black points correspond to (3 = 3.80, dark magenta points correspond to /3 = 3.90, green points correspond to /3 = 4.05 
and blue points correspond to (3 = 4.20 (see Appendix^). The dashed horizontal line correspond to the experimental value of 



AJ 2 



M 2 . The physical value of the symmetric combination of the light quark masses is rh u d{MS, 2 GeV) = 3.6(2) MeV. 



Using eq. (71), the pion mass difference extracted from the fits of Figure [I] can be converted in the 



physical result for — M^ multiplying the results by the factor 

e 2 = e 2 = 47TQW = (72) 
and by two powers of the inverse lattice spacing a determined within the isosymmetric theory (see Ta- 



ble |Aj). Indeed, as shown in the diagrammatic analysis of section 5.2, the pion mass difference is a genuine 
isospin breaking effect and the change of the lattice spacing as well as the change of the average up- 
down quark mass enter at higher orders in the perturbative expansion. In Figure [2] we show the data for 
M 2 + — M 2 converted in physical units at the four different values of lattice spacings used in this paper and 
for the different values of rh u d used in the simulations. The horizontal black line in the figure corresponds 
to the experimental determination of the pion mass difference squared and the lattice data are very close 
to the experimental value. 
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Our pions are heavier than the physical ones and our lattice data need to be extrapolated toward 
the chiral limit. Furthermore, QED is a long range interaction and we have to cope with the associated 
power-law finite volume effects. The chiral extrapolation and the removal of lattice artifacts from simulated 
numerical data is the subject of section [9] 



7 Tuning critical masses 



In order to extract physical informations from the expression for M K + — M K o, see eq. (69) above, we 
first need to obtain a numerical determination of the electromagnetic shift of the critical masses of the 
light quarks. Our results for the kaon mass splitting have been obtained within the electro-quenched 
approximation that, consistently, we employ in this section to calculate Am^. 

As discussed in section [4~T| we can use two different conditions to obtain a numerical estimate of the 
divergent parameters Am^. The first strategy, based on Dashen's theorem, consists in imposing the 
validity of the continuum SU (3) chiral limit relations 

lim M^o = lim M K o = , (73) 

m/i— >0 rh/i— >-0 

where rhf = {m n ,m^,m s }. Relying on the determination of the QCD critical mass performed in 
ref. [l4] and using eqs. ([65]) and (68), we have that in the electro-quenched approximation the neutral pion 



and neutral kaon masses vanish for AmJ r given by 




Araf = — e l lim — — , (74) 



where / = {u,d,s}. From the numerical point of view, the parameters Am^ r have to be determined as 
accurately as possible because they are needed in order to cancel a linear ultraviolet divergence present 



in eq. (69). The numerical problem with eq. (74) is that the associated determination of AmJ r requires a 
chiral extrapolation and this in turn introduces larger uncertainties compared to the alternative method 
discussed in section |4.1[ namely the numerical determination of the electromagnetic critical masses based 



on the use of the WTI of eq. (41) 



By applying the methods of section [5] to the Ward-Takahashi identity Wf(g) = 0, i.e. by applying the 
differential operator A to the full theory parity-odd correlator (l.h.s. of eq. ([41])) 



+ 



W f 



(£) = -V <(^> =-VoTr{ 1 ° S+[U,A;g;t,p = 0]^ Sj[U,A;g;-t,p = 0] }=0, (75) 
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Figure 3: Left panel: determination of Amy according to eq. |76| > for the simulation corresponding to /3 = 4.20 and 
(am u d)° = 0.0020 (see Appendix [a} . As expected the combination of correlators appearing into eq. \76\ give a constant 
plateau in time from which we extract AmJ r . Right panel: numerical results for AmJ r for the different simulations. Black 
points correspond to (3 = 3.80, dark magenta points correspond to (3 = 3.90, green points correspond to (5 = 4.05 and blue 
points correspond to (3 = 4.20. As expected the critical mass counter-terms depend very mildly from the simulated symmetric 
light quark mass (am u d) : the small dependence is due to statistical fluctuations and (small) cutoff effects. 



one obtains the following alternative definition of Am' 



AW f = 



Am , 



V 




V 



(76) 



Note that the two definitions of eqs. (74) and (76) have the same "structure" in terms of corrected cor- 



relators. Indeed, the Dashen's theorem is a consequence of the chiral WTI of the continuum theory and, 
concerning valence flavour doublets, eq. (41) is the chirally twisted version of one of these relations. From 



the numerical point of view, however, the great advantage of eq. (76) with respect to eq. (74) is that the 
first does not require chiral extrapolations. 



In the left panel of Figure|3]we show the combination of correlators appearing into eq. (76) as a function 
of time for the simulation at (3 = 4.20 and (am u d)° = 0.0020, see Appendix [aJ As expected, coming from 
a WTI, the numerical data exhibit a very long plateau from which we obtain a reliable determination 
of AmJ. We have similar results for the other values of quark masses and lattice spacings simulated in 
this paper. In the right panel of the same figure we show, for each lattice spacing, AmJ r as a function 
of (am u d)°. As expected the results are almost insensitive to m Q ud . The tiny dependence that can be 
appreciated in the figure is due to statistical fluctuations and to (small) cutoff effects. In the following we 
use for each simulated value of m Q ud the corresponding determination of AmJ r and subsequently extrapolate 
the results for 



M| 



to continuum, infinite volume and chiral limits, see section [9J 



3.90 where the number of simulated values of do allow 



As a cross-check of our results, at (3 
a reliable chiral extrapolation, we have compared the determination of AmJ r obtained by using eq. (74) 
with the values in Figure [3] Though the two determinations may differ because of cutoff effects, the two 
numerical values agree within the errors (that are two orders of magnitude larger for the value of AmJ r 
obtained from the chiral extrapolation). 
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8 Kaon masses and separation of QED from QCD LIB effects 



With a reliable numerical determination of Am^ r and Am^ r we are now in the position of using the kaon 



mass difference formula of eq. (69) for physical applications. By defining 



R 



K 



O' 



K 



o 



T>exch 

H K - 



<S> 

o 



r>self 
K 




O 



(77) 



in the electro-quenched approximation we have 



M K+ -M K o = -2Am ud d t R% - (Am c J - Amf) d t R k K + (e 2 u - ej)e z d t R<g ch - R s ° lf . (78) 



The kaon mass splitting is a physical quantity and, since the electric charge does not renormalize at 
first order in a 
M K+ 

in order to predict the mass splitting of other hadrons, as for example the neutron-proton mass difference. 



em , the right hand side of eq. (78) can be made finite and equal to the physical value of 
M K o by properly tuning the bare parameter Amy. Afterward, the parameter Aray can be used 



Eq. (78) can also be used for introducing a renormalization prescription to separate QED from QCD 
isospin breaking corrections. This separation is not needed when, as in this paper, simulations are per- 
formed in the full theory but it may be useful in practice, as discussed in our previous work on the subject 2 
or in ref. p]. To this end, we need to express eq. (78) in terms of the renormalized light quark masses. Note 
that the bare parameters m u d and Amy of the full theory mix under renormalization because the two 
light quarks have different electric charges and, consequently, different renormalization constants Z mu (/i) 
and Zm d (fjL). Specifically, we have 



Am ud 



m d 



Amy ray 



Z. 



ud 



(79) 



where we have defined 



1 

Zud 



2\Z, 



1 

%ud 



2\Z, 



(80) 



The mixing does not occur in the isosymmetric theory where the quarks are neutral with respect to 
electromagnetic interactions and we have 



1 

zL 



1 



(81) 
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In the maximally twisted mass regularization used in this paper = Zp, see refs. [9||T5[, and in our 
numerical analysis we have used the non-perturbative results for Zp((3°) obtained in ref."[17| and listed in 
Table |A| By neglecting all the contributions of 0(e 2 Amy), eq. (79) can be rewritten as 



Am ud = Z^ Am ud ■ 



m ud 



(82) 



The electromagnetic and strong contributions to the kaon mass difference are conveniently separated 
by plugging the previous expression into eq. (78) and by defining 

[M K+ - M K o] QED (fi) = -2m ud ^ - (Am- - AmJ)d t R k K + (e 2 u - e 2 )e 2 d t R% ch ~ R% lf 



[M K+ - M K of CD ( M ) = -2Am ud (Z^ d t R%) , 
M K+ - M K o = [M K+ - M K o] QED (fi) + [M K+ - M K o] QCD (fi) 
This prescription was used in ref. \2\ where we determined Am„d(/i). 



(83) 



Given a conventional separation between electromagnetic and strong isospin breaking corrections, vi- 
olations to Dashen's theorem are parametrized in terms of "small" parameters and, concerning the kaon 
mass difference, one has (see refs. [lj|2]) 



[il4 + - M% } QED ( M ) - [Ml + - Ml^ u ( M ) 
Ml + -Ml 



2 iQED 



e 7 (/x) 



[M 



2 

K+ 



M 2 K0 ] QED (») 



M 2 , - M 2 



0(a em Am ud ) 



(84) 



As observed in ref. [2] one can in principle specify a renormalization prescription by fixing a value for e 1 
and by using eq. ( [84] ) to compute the corresponding value of Z ud {\±). This is not the strategy followed 
in this paper. Here, by relying on the perturbative result in the MS scheme that can be extracted from 
ref. 



18 



Z ud (MS,fi) 



(e 2 d -ejy 
32tt 2 



6 log(a/i) - 22.596... 



7° 



(85) 



we employ eq. (84) to calculate e 7 (MS,/i). The estimate provided for this quantity in ref. [l] and used in 
ref. [2] to calculate the leading QCD isospin breaking corrections to the decay rate is e 1 ~ 0.7. 

In Figure |4] we plot the fits of the different ratios of corrected correlators defined into eqs. ( 77 ), performed 
according to eq. ( plj ). In Figure [5] we plot s 1 {MS,2 GeV) as obtained from our numerical simulations, 
by using our lattice results for both the numerator and the denominator of the ratio appearing in 
The data obtained at unphysical values of rh ud 
to confirm that e 

next section we discuss the chiral extrapolation and the removal of cutoff and finite volume effects from 
the data of Figure [5| 



i.e. 
eq. (184) 



at finite lattice spacing and at finite volume seem 
7 — 0.7 is a reasonable estimate for the Dashen's theorem breaking parameter. In the 
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Figure 4: Fits of K%(t) (top-left), R k K {t) (top-right), R e £ ch (t) (bottom-left) and R s ^ lf \t) (bottom-right) according 
to eq. \6l\ . The dark magenta points correspond to (3 = 3.90 and (am u d)° = 0.0030, the green points to /3 = 4.05 and 
(am u d) = 0.0060 while the blue points correspond to f3 = 4.20 and (am u d) = 0.0065 (see Appendix^). Data are in lattice 
units. 
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Fi gure 5: e 7 (M5, 2 GeV) as a function of rh u d as extracted from lattice correlators. The dashed horizontal line correspond 
to the value e 1 = 0.7 used into ref. 2 . Black points correspond to /3 = 3.80, dark magenta points correspond to (3 = 3.90, 
green points correspond to f3 = 4.05 and blue points correspond to /3 = 4.20 (see Appendix [A}. 
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9 Chiral extrapolations and lattice artifacts 



In this section we discuss the chiral extrapolation and the removal of cutoff and finite volume effects from 



M 2 



M 2 



(see section [6] and Figure |2| and from e 7 (see previous section and Figure [HJ. 



We rely on the chiral formulae of ref. 19 (see also ref. 5 ) where the authors studied the dependence 
of — M^o and — M^o on the quark masses together with finite volume corrections by using an 

effective chiral lagrangian at NLO with the inclusion of electromagnetic interactions. They defined QED 
on a finite volume by considering the same infrared regularization used in this paper. In the case of the 
pion mass difference the formulae corresponding to the rif = 2 theory with a quenched strange quark are 



f**[C,K\= [M 2 + -M 2 ] 



2e 2 F 2 



C - (3 + 4C) 



Ml 



32tt 2 F 2 



log ( M ) + A' (/i) 



f?\C\ = [M 2 + - M 2 ] (L) - [Ml + - M 2 ] (oo) 



4irL 2 



[F 2 (M T L)-4Cff 1 (M w L)] . (86) 



The corresponding formulae for £ 7 are 



•2ef 



3 
C 



Ml 
32tt 2 F 2 



Ml 
32tt 2 F 2 



log 



log 



M% 



Mi 



(87) 



and 



£ 7 (L) — £ 7 (oo) 



87rL 2 F 2 C 



H 2 (M K L) - H 2 {M V L) 

H\{MkL) - Hi(M w L) 



i 9 p sea | 9 sea 

8ttL 2 F 2 V 3 + " + d 



In the previous expressions = 2B^m u d and Mk = Bo(m s +rh u d) with .Bo and Fo the QCD low energy 
constants entering the leading order chiral lagrangian. Fq is normalized so that the decay constant of the 
physical pion is about 90 MeV. Furthermore, C is the single electromagnetic low energy constant entering 
the leading order lagrangian while K(/a), K 1 ^) and K 2 (ji) are combinations of electromagnetic low energy 
constants at next-to-leading order. Note that the formulae for e 1 depend upon the charges 
of the sea quarks that have to be set to zero in the electro-quenched approximation. 

The functions Hi^(x) entering the finite volume correction formulae are 



' a and e s d ea 



H 2 (x) = 



du — -rr^ e 



-2k 



I 



+ CO 



du 



S(u) 



x\fu erf ( x 



(89) 
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Fi gure 6; Dependence of the functions H\{x) and H2(x) with respect to x = kL. In our numerical simulations 
we have 3.2 < M^L < 5.8 and 5.1 < M^L < 7.0. The linear approximation of the function H2(x) is explicitly given by 
H 2 (x) ~ -k(2 + z). 



where 



F(x) = [03 (0, i/x)f - 1 , S(x) = x 3/2 - F(x) , k = / + du^- = 2.8373 . . . 

Jo u 



are expressed in terms of the Jacobi elliptic ^-function, 

# 3 (0, i/x) = ^3(0, ix) = I 1 + 2 e ~ 



(90) 



(91) 



The formulae defining the functions H\^{x) contain rather involved integral expressions but, in the 
range of values of x = M^,kL corresponding to the meson masses and lattice volumes simulated in this 
paper, i.e. 3.2 < x < 7.0, the function H2(x) is almost a linear function of x and H\{x) can be safely 
neglected (see Figure [6]). From the asymptotic expressions 



[M 2 + - M 2 ] (L) - [Ml + - Mlo] (oo) 
e 1 (L) — e 7 (oo) 



e 2 k ( M % , 2 \ 



2l - f — + 
4tt V L L 2 J 

k M K - M v 
87rF 2 C L 



(92) 



one recognizes that QED, a long range interaction, introduces sizable power-law finite volume effects on 
hadron masses that have to be eliminated by extrapolating numerical data obtained on different physical 
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Figure 7: Combined chiral, continuum and infinite volume extrapolations of M 2 



■AC 



The grey points are the data 



as extracted from lattice simulations and converted in physical units. The red point and the solid red curve is the result 
of the extrapolation. Black points correspond to /3 = 3.80, dark magenta points correspond to (3 = 3.90, green points 
correspond to j3 = 4.05 and blue points correspond to j3 = 4.20 (see Appendix [A}. The dashed horizontal line corresponds 
to the experimental value of — M 2 . For each panel the colored points have been obtained by subtracting the fitted 

lattice artifacts. Note the two points at rh u d ~ 22 MeV that have been obtained at (3 = 3.9 and at the same value of sea 
quark mass but on different volumes: within quoted errors the two points differ at finite volume (grey) and coincide after the 
removal of discretization effects (dark magenta). 



volumes. These are predicted to be as large as 30% in our simulations and to be larger at heavier meson 
masses. 

In order to extrapolate lattice data for the pion mass difference — M^ we have performed three 
fits differing in both the dependence on the average light quark mass m u d and the parametrization of the 
lattice artifacts. We have considered the following functions 



mC,K,A^B„] 



fx*[C,K] + fr[C} + A„ [a ]' 



C + Km ud +^-+A v [a ] 2 , 
C + K m ud + ^+A„ [a ] 2 . 



(93) 



The function /f correspond to the SU (2) chiral fit and has been obtained by using the chiral perturbation 
theory results of eqs. (86) for the dependence on rh u d and the finite volume corrections. The other 
two functions correspond to phenomenological fits performed assuming a linear rh u d dependence and 
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parametrizing finite volume effects with a term either proportional to 1/L or to 1/L 2 . In all cases we have 
added a term proportional to (a ) 2 in order to estimate cutoff effects. 

In Figure [7] we show the combined chiral, continuum and infinite volume extrapolations corresponding 
to the three fitting functions defined above. In each plot the solid red curve represents the fitted function ff 
evaluated at a = and at 1/L = 0, i.e. to the result of the continuum and infinite volume extrapolations. 
The grey points are our lattice results for M 2 + — M 2 , already shown in Figure [2] The red point is the 
result of the extrapolation at the physical value rh u d(MS,2 GeV) = 3.6(2) MeV determined within the 
isosymmetric theory in ref. [20] . The remaining colored points have been obtained from the corresponding 
grey points by subtracting out the lattice artifacts as determined by the fit. As the colored points are 
the results of the continuum and infinite volume extrapolations they have larger errors compared to the 
corresponding grey points. The three fits give consistent results for M 2 + — M 2 , though with different 
errors. In particular, the fit that include a 1/L finite volume term gives larger errors on the extrapolated 
points than the other fitting functions. By comparing the /f panel with the other two, we see that the 
finite volume effects obtained from the fits f£ and f£ are considerably smaller than the one-loop chiral 
perturbation theory prediction. Similar results for finite volume effects, parametrized by a 1/L 2 term and 
fitted to lattice data, have been obtained by the authors of ref. 4 . The value of x 2 /dof corresponding to 
the fit fl is 1.2, while for the fits /J and fl we have x 2 /dof = 1.0. 

We obtain our final estimate of the pion mass difference by taking as central value and as statistical 
error the results of the SU{2) chiral fit, i.e. /f , while we estimate our systematic error by taking half of 
the difference between the f£ and fg results. We get 

M 2 + - M 2 = 1.44(13)(16) x 10 3 MeV 2 . (94) 
This result compares nicely with the experimental determination 

[M 2 + - M 2 ] exp = 1.2612(1) x 10 3 MeV 2 , (95) 
suggesting, a posteriori, that the effect of having neglected the disconnected contribution of 0(a em m U( i) 



appearing in eq. (66) is smaller or of the same order of magnitude as the other uncertainties affecting our 
result. 

We now discuss the extrapolations of the results for e 1 {MS 1 2 GeV). In general we expect reduced 
lattice artifacts for e 1 with respect to M 2 + — M 2 because of possible cancellations of systematics effects 
between the numerator and the denominator in the ratio of eq. ( |84| . Within the quoted errors, the lattice 
data shown in Figure [5] are fairly flat in m u d so that we have not attempted a SU (3) chiral extrapolation 



using eqs. (87) and (88). Linear chiral extrapolations lead to vanishing slopes and errors on the fitted 
results of the same order of magnitude of the ones shown in Figure [8j corresponding to constant chiral 
extrapolations. More precisely, the fit functions shown in Figure [8] are 

fl[E,A £ ] = E + A £ [a } 2 , 
f![E,A £ ,B £ } = E+^+A £ [a } 2 , 

f§[E, A £ , B e ] = E+^+A £ [a ] 2 , (96) 
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Figure 8: Combined chiral, continuum and infinite volume extrapolations of e 1 (MS,2 GeV). The grey points are the 
data as extracted from lattice simulations. The red point and the solid red curve is the result of the extrapolation. Black 
points correspond to ft = 3.80, dark magenta points correspond to j3 = 3.90, green points correspond to ft = 4.05 and blue 
points correspond to f3 = 4.20 (see Appendix |Aj. For each panel the colored points have been obtained by subtracting the 
fitted lattice artifacts. 



In order to obtain an estimate of the systematic errors associated with our result, we have parametrized 
cutoff effects with a term proportional to (a ) 2 and finite volumes effects with terms vanishing as 1/L and 
as 1/L 2 . 

For each plot in Figure [SJ the solid red line corresponds to the fitted function ff evaluated at a = and 
at 1/L = 0, i.e. to the result of the continuum and infinite volume extrapolations. The grey points are our 
lattice results for s 1 {MS, 2 GeV) already shown in FigureJHJ The red point is the result of the extrapolation 
at the physical value m u d(MS,2 GeV) = 3.6(2) MeV determined in the isosymmetric theory in ref. [20]. 
The remaining colored points have been obtained from the corresponding grey points by subtracting out 
the lattice artifacts as determined by the fit. As the colored points are the results of the continuum and 
infinite volume extrapolations they have larger errors compared to the corresponding grey points. The 
three fits give consistent results for s 1 {MS, 2 GeV). The \ 2 /dof of all the three fits is 1.0. 

We obtain our final result of e 1 (MS, 2 GeV) by taking as central value the average of the maximum 
(/f ) an d minimum (/f ) central values and as statistical error the one obtained from the /| fit. In order to 
estimate the systematic error associated with the fits we take half of the difference between the /f and f± 
results (13%). By taking the deviation of our result for M 2 + — M 2 from the experimental value we obtain a 
rough estimate of the error associated with the neglected contributions, i.e. the M^o disconnected diagram 
and all the terms of 0(a em A , m u d) or higher, and add in quadrature a 15% uncertainty to the systematic 
error. The uncertainty associated with the electro-quenched approximation is estimated by using the 



chiral formula for e 7 given in eq. (87). More precisely, by taking C from the fit /f and by neglecting the 
variation of the terms K 1 and K 2 , we evaluate the ratio e 1 (e s ^ a + e s d ea = l/3)/£ 7 (e^ ea + e s d ea = 0) and add 
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in quadrature the resulting 12% uncertainty to the systematic error. We get 



£ 7 (M#, 2 GeV) = 0.79(18)(18) . (97) 

This result corresponds to the following separation of QED from QCD isospin breaking corrections to the 
kaons mass difference (see eqs. ([83]) above) 



[M%+ - M 2 K o] QED (MS, 2 GeV) = 2.26(23)(23) x 10 3 MeV 2 , 

[M 2 K + - M 2 KQ ] QCD (MS, 2 GeV) = -6.16(23)(23) x 10 3 MeV 2 , (98) 

where the QCD contribution is obtained by imposing the experimental constraint 

[M 2 K+ - M 2 KQ ] exp = -3.903(3) x 10 3 MeV 2 . (99) 

The experimental input for the kaon mass splitting also allows a determination of the up-down mass 



difference and, from the second of eqs. (83), we obtain 



[m d - m u ](MS, 2 GeV) = 2Am ud (MS, 2 GeV) = 2.39(8) (17) MeV , 

rri 

-^(MS,2GeV) = 0.50(2)(3), (100) 

m d 

The results for the light quark mass ratio has been obtained by combining the determination of Am u d 
performed in this paper with the result m u d(MS, 2 GeV) = 3.6(2) MeV obtained in ref. [201. The authors 
of ref. 20 have extracted the symmetric light quark mass the strange quark mass and the 



lattice spacing a by performing simulations of the isosymmetric theory and using the necessary number 
of hadronic inputs to calibrate the lattice. By following this procedure their result for rh? ud differs from 
rhud-, as defined in eqs. (Js]) at /i = 2 GeV, by isosymmetric 0((e^ + e 2 d )a ern ) contributions. These terms 
do not affect our results for the quark mass ratios since m u /rhd is a function of Arh u d/rh u d only and, in 
turn, 

Am ud Am ud , A A , , . 
~ — o 1- 0{a em Am ud ) , (101) 



m ud m 



ud 



but represent a systematic error for the determinations of rh u and rhd which follow from eqs. (100), namely 

m u (MS, 2 GeV) = 2.40(15)(17) MeV , 

m d (MS, 2 GeV) = 4.80(15)(17) MeV . (102) 
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The argument used in the case of rh u /rhd also applies to the flavour symmetry breaking parameters 



R(MS, 2 GeV) 



Q(MS,2 GeV) 



~rh s 


- rhud 


_ rhd 


- rhu 


\rh 2 s - 


~Kd 


™d 


-ml 



(MS, 2 GeV) = 38(2) (3) , 



(MS, 2 GeV) 



23(1)(1) , 



(103) 



that have been calculated starting from the relations 



_1_ _ 2Am ud 
R rhP Q — rhP ir 



0(a ern Arh ud ) 



1_ _ 4Amy rh Q ud 
r ~ (m2) 2 -K,) 2 



0(a ern Arh ud ) 



(104) 



with the result rh° s (MS, 2 GeV) = 95(6) MeV obtained in ref. [20]. 

Using our own determination of e 1 we are here in position of updating the results obtained in ref. [2] 
for the QCD isospin breaking effects on the decay rate and the neutron-proton mass difference. We 
get 



1 QCD 



1 



(MS, 2 GeV) 



-0.0040(3)(2) , 



Fk/F^ 

[M n - M P ] QCD (MS, 2 GeV) = 2.9(6)(2) MeV , 



(105) 



where F K + and F n + are the charged kaon and charged pion decay constants in QCD with rhd 7^ and 
Fk and F n are the corresponding quantities in the isosymmetric theory. Analogously, M n and M p are 
respectively the masses of the neutron and of the proton in QCD with rhd ^ ™>u- The results in eqs. (105) 
have now been obtained from a self consistent non-perturbative lattice calculation. 



10 Conclusions 



In this paper we have shown that leading isospin breaking effects on hadron masses can be conveniently 
calculated on the lattice by starting from simulations of the underlying isosymmetric QCD theory and by 
expanding the path-integral in powers of a em and rhd — ™u- We have discussed all the details necessary 
for applying our method to the calculation of isospin breaking corrections to any observable and discussed 
the renormalization of the corrected correlation functions. 

In particular, we have shown how the ultraviolet divergences generated by the contact-terms of the two 
electromagnetic currents can be absorbed in a redefinition of the bare parameters of the full theory with 
respect to the corresponding values of isosymmetric QCD. We have also shown that the linear divergences 
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associated with the shift of the critical masses of the quarks, a problem to be addressed when the fermion 
action is discretized by using a Wilson term, can be subtracted out by determining the associated counter- 
terms with high numerical precision. 

By using the proposed method we have derived theoretical predictions for the pion mass splitting and 
for the up and down quark masses. We have also implemented a well defined renormalization prescription 
to separate QED from QCD isospin breaking corrections to hadron masses allowing to determine the 
electromagnetic contribution to the kaon mass splitting and the associated value of the Dashen's theorem 
breaking parameter e 1 . These results have been used in order to update our previous determinations [2] 
of the QCD contributions to the neutron-proton mass splitting and the decay rate. 

The results obtained in this paper are affected by systematic errors. Particularly important are those 
associated with the chiral extrapolation required because our pions are heavier than the physical ones. 
Another important source of systematics errors comes from finite volume effects. These are not peculiar to 
our method. We estimate finite volume effects by using the results of effective field theory calculations and 
by fitting them to lattice data obtained at different physical volumes. The finite volume effects arising from 
the fits turn out to be considerably smaller than chiral perturbation theory predictions, though we cannot 
make this statement more quantitative until we have results at the physical pion masses and on larger 
physical volumes. Finally, our results have been obtained in the rif = 2 theory and neglecting certain quark 
disconnected diagrams. Concerning the pion mass splitting, we have neglected the disconnected diagram 



appearing in eq. (66), an 0(a era rh u d) correction to M n o that phenomenologically is expected to be of the 
same order of magnitude of the other 0(a em Am u d) contributions neglected in this paper. The results 
for the kaon mass difference have been obtained by relying on the electro-quenched approximation that 
consists in treating dynamical quarks as electrically neutral particles. We have provided all the formulae 
necessary to remove these approximations. This will be the subject of future work. 
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Appendix A Gauge ensembles 



In this work we have used the n f = 2 dynamical gauge ensembles generated and made publicly available by 
the European Twisted Mass Collaboration (see Table |A|) . These gauge configurations have been generated 
within the isosymmetric theory, by using the so-called Twisted Mass lattice discretization of the QCD 
action [9|fl0] , see eqs. (|30|) and (32). For the different gauge ensembles used in this work, the values of 

(ref. [l4]), lattice spacing a (ref. [20]), strange valence 
) are given in Table [aJ 



the critical hopping parameter fco = l/(2raQ r -\ 
quark mass (am s )° (ref. [20]), renormalization constant Z P (ref. 



17 



In the following we give a dictionary to translate in the operator language the diagrammatic notation 
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k 



(am ud )° (am s )° L/a N conf q° (fm) Z° P (MS,2 GeV) 



3.80 0.164111 



0.0080 0.0194 24 240 0.0977(31) 0.411(12) 

0.0110 24 240 



3.90 0.160856 



0.0030 
0.0040 
0.0040 
0.0064 
0.0085 
0.0100 



0.0177 



32 
32 
24 
24 
24 
24 



150 
150 
240 
240 
240 
240 



0.0847(23) 0.437(07) 



4.05 0.157010 



0.0030 
0.0060 
0.0080 



0.0154 



32 
32 
32 



150 
150 
150 



0.0671(16) 0.477(06) 



4.20 0.154073 



0.0020 
0.0065 



0.0129 



48 
32 



100 
150 



0.0536(12) 0.501(20) 



Table A: Gauge ensembles used in this work. The gauge configurations have been generated within the isosymmetric theory 
with Nf = 2 dynamical flavours of maximally twisted quarks of mass (am u d)° . The strange quark mass (am s )° has been 
used for valence propagators. The hopping parameter ko is related to the critical mass parameter mg r appearing in the main 
body of the paper by the relation ko = l/(2mg r + 8). 

used in the main body of the paper. To this end, it is convenient to define the following local operators 



v ++ 



(x) 



fg 



(x) = i 



T tg + 



(x) = j,+ {x) 7m U,{x)^{x + ^) + ^(x + m ) +n5 + ^ Ul(x)^(x) , 



l fg 
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(106) 



In terms of them we can now translate in local operator language the diagrammatic relations appearing 
in the main text. For example 



+ 



(107) 



where, as a general rule, if two different quark lines have the same color they have the same mass (i.e. 
mi = 7712). As further examples we have 



<^> = £T<P^)S 2 - 3 -(2/)P 3 -i + (0)}, 



<^> = i^T(P+-( a; )P 2 -3-(y)P 3 - 1 + (0)). 



(108) 



Examples involving the insertion of the electromagnetic currents are given here below 




£T(P+-(z) [V^iy) [V z - A -] v {z)P^{0))D^z), 



-<^> = £T(P+-(z) [V 2 3-]"(y)i£ + (0) [V+ + ]"(z))D ltv (y,z) 




±J2 T ( P i + 2 "W ^23 + (0) [T 3 + + f (y) ) D^(y,y) . 



(109) 



We finish with an example concerning a quark disconnected diagram, 




^T(P+-(x)P 2 -3 + (0) [V+ + ]"(y) [Vt i ± Y(z))D flv (y,z). (110) 



yz 



All the other diagrams appearing in the text can be translated into the operator language following the 
correspondence illustrated above. 
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